rm(list=ls())

# Load packages
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(readxl)
library(ggplot2)
library(gridExtra)
library(readstata13)
library(haven)
library(dplyr)
library(forcats)
library(ggeasy)
library(gt)
library(tidyverse)
library(broom)
library(RDHonest)

# Load placebo data
cep_all_placebo <- read_dta("Data/cep_all_placebo.dta")
dosage_ranges <- c(1, 2, 3, 4, 5)
sample_list <- lapply(dosage_ranges, function(i) subset(cep_all_placebo, treatment_dosage_coup >= -i & treatment_dosage_coup <= i))
names(sample_list) <- paste0("sample", dosage_ranges)

# Clustered SE function
cluster_se <- function(model, data, type = "HC1") {
  idx <- as.integer(rownames(model.frame(model)))
  vc  <- sandwich::vcovCL(model, cluster = data$cluster[idx], type = type)
  lmtest::coeftest(model, vcov. = vc)
}

# Model runner
run_models <- function(samples, formula, dv_name, table_title, dep_caption, cov_label) {
  models_lm <- lapply(samples, function(s) lm(formula, data = s))
  models_cl <- Map(cluster_se, models_lm, samples)
  sample_sizes <- sapply(models_lm, nobs)
  r_squared <- sapply(models_lm, function(m) summary(m)$r.squared)
  mean_dv <- sapply(models_lm, function(m) mean(m$model[[dv_name]], na.rm = TRUE))
  
  stargazer(models_cl,
            type = "latex",
            title = table_title,
            keep = "^treat1$",
            column.labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
            covariate.labels = cov_label,
            dep.var.caption = dep_caption,
            dep.var.labels = dv_name,
            dep.var.labels.include = FALSE,
            add.lines = list(
              c("Obs.", sample_sizes),
              c("R^2", round(r_squared, 3)),
              c("Mean DV", round(mean_dv, 3))
            ),
            align = TRUE)
}

# Table D1: Left Identification, Less than High School
tabd1 <- run_models(sample_list,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table D1: RD Estimates 1973 Coup on Left-wing Ideological Identification (Sample: Less than High School)",
           dep_caption = "Dependent Variable: Left Identification",
           cov_label = "Coup")

print('Table D1 complete')

# Table D2: Hard Left Identification, Less than High School
tabd2 <- run_models(sample_list,
           left_hard ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left_hard",
           table_title = "Table D2: RD Estimates 1973 Coup on Hard Left-wing Ideological Identification (Sample: Less than High School)",
           dep_caption = "Dependent Variable: Hard Left Identification",
           cov_label = "Coup")

print('Table D2 complete')

# Load data with controls
cep_all <- read_dta("Data/cep_all.dta")
sample_list <- lapply(dosage_ranges, function(i) subset(cep_all, treatment_dosage_coup >= -i & treatment_dosage_coup <= i))
names(sample_list) <- paste0("sample", dosage_ranges)

# Table D3: Left Identification w/ Controls
tabd3 <- run_models(sample_list,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta) + sexo + edad,
           dv_name = "left",
           table_title = "Table D3: RD Estimates 1973 Coup on Left-wing Ideological Identification With Control Variables",
           dep_caption = "Dependent Variable: Left Identification",
           cov_label = "Coup")

print('Table D3 complete')

# Table D4: Hard Left Identification w/ Controls
tabd4_hardleft <- run_models(sample_list,
           left_hard ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta) + sexo + edad,
           dv_name = "left_hard",
           table_title = "Table D4: RD Estimates 1973 Coup on Hard Left-wing Ideological Identification With Control Variables",
           dep_caption = "Dependent Variable: Hard Left Identification",
           cov_label = "Coup")

# Table D4 (continued): Soft Left Identification w/ Controls
tabd4_softleft <- run_models(sample_list,
           left_soft ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta) + sexo + edad,
           dv_name = "left_soft",
           table_title = "Table D4: RD Estimates 1973 Coup on Soft Left-wing Ideological Identification With Control Variables",
           dep_caption = "Dependent Variable: Soft Left Identification",
           cov_label = "Coup")

print("Table D4 complete")

# Load alternative dataset (Bicentenario)
bic_all <- read_dta("Data/bic_all.dta")
sample_list <- lapply(dosage_ranges, function(i) subset(bic_all, treatment_dosage_coup >= -i & treatment_dosage_coup <= i))
names(sample_list) <- paste0("sample", dosage_ranges)

# Table D5a: Left Identification, Bicentenario
tabd5_left <- run_models(sample_list,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table D5: RD Estimates 1973 Coup on Left-wing Ideological Identification (Alternative Data)",
           dep_caption = "Dependent Variable: Left Identification",
           cov_label = "Coup")

# Table D5b: Right Identification, Bicentenario
tabd5_right <- run_models(sample_list,
           right ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "right",
           table_title = "Table D5: RD Estimates 1973 Coup on Right-wing Ideological Identification (Alternative Data)",
           dep_caption = "Dependent Variable: Right Identification",
           cov_label = "Coup")

# Table D5c: Center Identification, Bicentenario
tabd5_center <- run_models(sample_list,
           center ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "center",
           table_title = "Table D5: RD Estimates 1973 Coup on Center Ideological Identification (Alternative Data)",
           dep_caption = "Dependent Variable: Center Identification",
           cov_label = "Coup")

# Table D5d: No Ideology, Bicentenario
tabd5_none <- run_models(sample_list,
           no_ideology ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "no_ideology",
           table_title = "Table D5: RD Estimates 1973 Coup on No Ideological Identification (Alternative Data)",
           dep_caption = "Dependent Variable: No Ideology",
           cov_label = "Coup")
